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Reciprocal  Distance  Squared  Method 
A  Computer  Technique  for  Estimating  Areal  Precipitation' 


Tsong  Chang  Wei  and  J.  L.  McGuinness2 

INTRODUCTION 

Determination  of  watershed  rainfall  is  basic  to  the  calculation  of  other  such  hydrologic  characteristics  as  runoff  and 
infiltration.  The  accuracy  of  areal  rainfall  estimation  inevitably  influences  the  result  of  these  other  computations.  This 
report  presents  a  computer-oriented  technique,  Reciprocal  Distance  Squared  (RDS),  developed  for  the  evaluation  of 
areal  precipitation.  The  method  imposes  a  grid  system  on  the  area  and  calculates  the  amount  of  rainfall  at  each  grid 
point  by  summing  the  products  of  gage  weight  and  measured  rainfal  at  nearby  rain  gages.  Gage  weights  are  computed 
as  fractions  of  the  reciprocal  distance  squared  between  the  point  and  the  rain  gages.  The  amount  of  areal  precipitation 
can  be  evaluated  by  integrating  the  grid  point  rainfall  within  the  boundary  of  the  area.3 

The  methods  commonly  used  in  calculating  areal  rainfall  are  the  arithmetic  mean,  Thiessen,  and  the  isohyetal 
method.  Among  these,  the  isohyetal  method  is  generally  considered  the  best,  and  the  accuracy  of  the  RDS  method  is 
judged  on  its  ability  to  reproduce  isohyetal  results.  Since  the  true  areal  rainfall  is  unknown,  no  objective  method  can 
be  adopted  for  determining  the  superiority  of  the  methods. 

To  reduce  the  laborious  manual  work  and  provide  reproducible  results,  computer  programs  for  calculating  Thiessen 
weights  and  for  drawing  isohyetal  maps  have  been  presented  by  Diskin  (3,4)4  and  Diskin  and  Davis  (5).  Diskin  (3) 
first  suggested  a  Monte  Carlo  method  for  evaluating  Thiessen  weights.  Later  (4)  he  proposed  a  rectangular  grid  system 
method  to  improve  the  result.  In  the  program  to  draw  isohyetal  maps,  Diskin  and  Davis  (5)  constructed  a  network  of 
lines  connecting  adjacent  rain  gages  and  then  used  linear  interpolation  to  find  coordinates  of  isohyetal  points  between 
two  rain  gage  stations.  By  connecting  points  with  equal  amounts  manually,  the  isohyetal  map  was  drawn. 

Rapid  progress  in  the  application  of  watershed  modeling  to  solve  hydrologic  and  water  resources  problems  has 
created  a  demand  for  an  effective,  precise,  and  computer  oriented  means  of  evaluating  areal  rainfall.  Several  new 
approaches  have  been  presented  in  recent  years  (1,  7,  10,  11).  The  general  procedures  of  these  approaches  can  be 
defined  as  follow:  A  rainfall  distribution  equation  is  derived  from  the  location  of  rain  gages  and  rainfall 
measurements.  A  grid  system  is  imposed  on  the  watershed  and  rainfall  at  all  grid  points  is  calculated  from  the 
equation.  Integrating  rainfall  at  grid  points  with  respect  to  the  area  within  the  watershed,  divided  by  the  watershed 
area,  gives  the  average  areal  rainfall. 

Utilizing  a  grid  system  to  supplement  hydrologic  data  where  no  measurement  has  been  conducted  is  known  as  the 
grid  square  technique,  and  it  has  been  used  by  Solomon  and  others  (11 )  to  estimate  precipitation,  temperature  and 
runoff  at  Newfoundland,  Canada.  Later,  Pentland  and  Cuthbert  (10)  applied  the  technique  in  generating  synthetic 
steamflow  of  ungaged  streams  in  the  New  Brunswick  region  of  Canada.  One  strong  advantage  of  the  technique  is  the 
capability  of  using  digital  computers  which  can  handle  huge  amounts  of  data.  Most  of  the  calculations  for  grid  points 
can  be  performed  as  a  matrix  operation  without  any  difficulty  in  a  computer. 


1  North  Appalachian  Experimental  Watershed,  North  Central  Region,  Agricultural  Research  Service,  U.S.  Department  of 
Agriculture,  Coshocton,  Ohio  in  cooperation  with  the  Ohio  Agricultural  Research  and  Development  Center,  Wooster,  Ohio. 

2  Formerly  Agricultural  Engineer  and  Research  Statistician. 

3  A  similar  technique  is  described  by  W.  T.  Sittner  in  an  unpublished  manuscript  and  outlined  in  "National  Weather  Service  River 
Forecast  System-Forecast  Procedures,"  NOAA  Tech.  Memo.  NWS  HYDRO-14,  1967. 

4  Italic  numbers  in  parenthesis  refer  to  Literature  Cited,  p.  13. 
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Many  schemes  have  been  suggested  for  determining  the  pattern  of  rainfall  on  a  watershed  area.  After  reviewing 
area-depth  formulas,  Court  (2)  suggested  the  bivariate  Gaussian  (or  normal)  distribution  which  provides  elliptical 
isohyets  as  a  means  of  describing  the  rainfall  distribution  pattern  of  storms.  Assuming  a  polynomial  distribution, 
Mandeville  and  Rodda  (7)  suggested  the  use  of  a  trend-surface  analysis  to  determine  the  best  fit  pattern.  When  the 
polynomial  form  was  used  in  the  distribution  equation,  Chidley  and  Keys  (1)  showed  with  a  rectangular  shaped 
watershed  that  the  Thiessen  type  gage  weight  can  be  used  to  simplify  the  computation  processes  of  areal  rainfall. 

There  are  some  shortcomings  in  the  approaches  just  outlined.  First,  the  distribution  pattern  of  the  rainstrom  is 
expressed  in  a  mathematical  form  that  could  never  be  obtained  precisely  in  actual  storms.  At  best,  the  distribution 
equation  indicates  the  averaged  tendency  of  rainfall  distribution  of  a  storm  but  not  the  real  distribution  pattern. 
Second,  it  is  necessary  to  determine  the  distribution  equation  for  every  rainstorm.  Third,  the  calculation  processes  are 
irreversible,  i.e.,  when  a  distribution  equation  has  been  established,  it  will  usually  not  estimate  the  measured  rainfall 
exactly  at  any  of  the  rain  gages.  Finally,  the  distribution  equation  fits  the  overall  trend  of  the  rainfall  pattern  and 
smooths  out  local  variations  that  may  be  important  in  the  input  of  a  watershed  model. 

Taking  into  account  the  deficiencies  mentioned  above,  the  RDS  method  considers  the  amount  of  rainfall  at  any 
ungaged  point  as  a  function  of  measured  rainfall  and  distance  from  the  point  to  the  nearby  rain  gages  alone.  Thus,  no 
rainfall  distribution  equation  is  required  in  the  calculation.  By  limiting  the  affecting  gages  to  the  closest  few  gages  and 
by  assigning  greater  weight  to  the  nearest  gage,  this  method  minimizes  the  tendency  to  smooth  out  the  rainfall 
distribution  pattern. 

Two  computer  programs,  written  in  Fortran  language,  were  used  for  the  computations.  The  examples  were  run  on  a 
terminal,  at  the  North  Appalachian  Experimental  Watershed,  Coshocton,  Ohio,  which  had  access  to  General  Electric's 
former  Mark  II  computer  time-sharing  service.5  Results  are  compared  to  those  calculated  from  conventional  methods. 
It  is  also  shown  that  with  a  proper  arrangement  of  the  proposed  method,  Thiessen  type  gage  weight  can  be  calculated 
from  the  geometric  factors  of  gage  location  and  from  the  grid  system. 

RECIPROCAL  DISTANCE  SQUARED  (RDS)  METHOD 

If  a  number  of  rain  gages  exist  on  a  watershed  and  an  estimate  of  the  amount  of  rainfall  at  an  ungaged  point  in  the 
area  is  required,  it  is  rational  to  make  the  estimate  using  the  gages  nearest  the  ungaged  point  and  neglecting  those  at 
remote  distances.  Even  among  nearby  gages,  different  factors  such  as  the  relative  distances  between  gages  and  the 
point,  altitude  of  each  gage,  and  the  slope  and  orientation  of  topography  around  each  gage,  results  in  each  gage  having 
a  different  true  weight  in  the  estimation.  But  in  most  cases,  the  strongest  weight  is  from  the  nearest  gage  and  weights 
reduce  proportionately  as  the  distance  increases.  The  following  assumptions  were  made  to  estimate  the  rainfall  at  an 
ungaged  point. 

First,  impose  a  Cartesian  coordinate  grid  system  on  the  watershed  so  that  the  entire  watershed  area  lies  in  the  first 
quadrant.  Then  assume  that  the  rainfall  at  any  ungaged  point  (x,  y)  in  a  watershed  can  be  determined  by  n  number  of 
nearby  rain  gages  around  the  point  where  n  is  less  than  the  number  of  rain  gages  available  in  the  watershed.  Also, 
assume  that  the  rainfall  at  this  point  is  porportional  to  the  rainfall  measured  at  the  n  rain  gages  and  inversely 
proportional  to  the  distances  between  the  point  and  the  rain  gages.  After  testing  the  reciprocal  distance  in  linear, 
square  and  cubic  forms,  the  isohyetal  map  of  the  square  form  was  found  to  most  closely  resemble  that  drawn  by 
conventional  methods.  Thus  using  the  squared  form,  the  amount  of  rainfall,  R,  at  point  (x,  y)  can  be  expressed  by 

Ft  =  &  Pj/Dj2)  /d  1/Dj2)  (i) 
i  =  1  i  =  1 

where  Pj  is  the  measured  rainfall  at  gage  i,  D{  is  the  distance  between  point  (x,  y)  and  gage  i,  and  n  is  the  number  of 
rain  gages  used  in  determining  rainfall  at  the  point  (x,  y). 

If  the  point  is  close  or  coincides  with  one  of  the  gages,  say  gage  j,  the  reciprocal  distance  squared  term  of  the  gage 
j  in  Equation  1  becomes  so  large  that  other  terms  in  the  summation  are  trivial  and  negligible.  Then,  since  the 
denominator  and  numerator  of  Equation  1  each  contains  the  same  term  of  the  reciprocal  distance  squared,  they 
cancel  one  another,  and  rainfall  at  point  R  becomes  identical  to  that  of  gage  j. 


sUse  of  a  trade  name  does  not  imply  endorsement  by  the  U.S.  Department  of  Agriculture.  It  is  used  only  for  the  purpose  of 
providing  specific  information  to  readers. 
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The  watershed  boundary  can  be  considered  as  a  series  of  straight  line  segments.  Then  the  area  of  a  watershed  can  be 
calculated  by  summing  up  numerous  small  polygons  of  area  AA.  if  the  apexes  of  successive  angles  in  a  polygon  are 
assigned  numbers  from  1  to  n,  the  coordinates  of  the  apexes,  (x,  ,  y,  ),  (x2 ,  y2  ),  - — ,  (xn,  yn)  can  be  read  from  the  grid 
system.  The  polygonal  area,  AA,  can  be  calculated  from 

n-1 

AA=  (x1Vn -xnY1>  +  2    (xj+1  yj-xj  yj+1)  (2) 
i=1 

where  x  and  y  are  the  x  coordinate  and  y  coordinate  of  the  apexes,  respectively,  and  the  subscripts  denote  the 
corresponding  number  assigned  to  each  apex. 

From  Equation  1,  the  depth  of  rainfall  at  all  points  of  each  polygon  can  be  calculated.  If  the  area  AA  is  chosen  small 
enough,  it  can  be  assumed  that  the  average  rainfall  in  the  area  AA  is  a  simple  arithmetic  mean  of  the  rainfall  at  the 
apexes  of  the  polygon.  That  is 

m 

AP  =  2  Rj/m  (3) 
j=1 

where  AP  is  the  average  depth  of  rainfall  in  an  area  AA,  Rj  is  the  calculated  depth  of  rainfall  at  the  apex  j,  and  m  is  the 
number  of  apexes  that  compose  the  polygon.  The  volume  of  rainfall,  AV,  in  the  area,  AA,  is  given  by 

m 

AV  =  aPaA  =  (2  Rj/m)  AA  (4) 
j=l 

if  AP  is  replaced  by  Equation  3. 

If  the  volume  of  rainfall  on  a  watershed,  V,  and  the  area  of  the  watershed,  A,  are  given,  the  average  depth  of  area 
rainfall,  P,  on  the  watershed  is 

P=V/A.  (5) 

Since  the  volume  of  rainfall  in  an  area,  AA,  is  given  in  Equation  4,  the  volume  of  rainfall,  V,  on  the  watershed  can  be 
obtained  by  a  summation  of  all  volumes  in  the  meshes.  Therefore,  Equation  5  can  be  rewritten  as 

k 

P  =  2  AVp/A  (6) 
p=1 

or 

k  mp 

P  =  [  2  (2    Rj/rrip) AA  p]  /A  (7) 
p=1  j=1 

where  k  is  the  number  of  mesh  areas  in  the  watershed  and  mp  and  AAp  are  the  number  of  apexes  and  the  area  of 
mesh  p,  respectively. 

The  amount  of  rainfall,  Rj,  in  Equation  7  at  each  apex,  j,  is  calculated  by  Equation  1.  Substituting  in  Equation  7, 

k     171  p        n  n 
P=  2  <2    [  (  2Pi/D|2)  /  (  2  1/Dj2)  ]  j/mp>AAp/A.  (8) 
P=1j=1       i=1  i=1 

If  Equation  8  is  completely  expanded,  it  can  be  seen  that  every  term  consists  of  P^  multiplied  by  a  function  of  Dj, 
mp,  AAp  and  A.  Considering  any  area,  AAp,  in  the  watershed,  the  number  of  corner  points,  mp,  is  a  fixed  number. 
Also  the  distances  from  the  corner  points  to  any  of  the  rain  gages  are  constants.  The  sum  of  reciprocal  of  distance 

n 

squared  in  the  innermost  parenthesis,  2    1/Dj2 ,  is  a  constant.  The  coefficients  of  all  P[  terms,  therefore,  become 

i=l 

constants.  Combining  similar  Pj  terms,  Equation  8  can  be  rewritten  as 

P=  2  WiPj  (9) 

i=l  v  ' 

where  Wj  is  a  function  of  the  geometric  factors,  Dj,  mp,  AAp  and  A,  which  are  all  constants.  Equation  9  actually  is  of 
the  same  form  as  the  Thiessen  weight  method,  i.e.  the  area!  rainfall  is  the  summation  of  the  products  of  gage  weight 
and  measured  rainfall.  Once  the  weights  of  the  rain  gages  have  been  calculated,  the  areal  rainfall  can  be  obtained  easily 
even  with  a  desk  calculator. 
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DESCRIPTION  OF  THE  COMPUTER  PROGRAMS 


Because  of  limitations  on  the  capacity  of  the  available  computer,  two  programs  were  written.  The  first  program 
(Program  A)  provides  grid  point  rainfall  for  drawing  an  isohyetal  map,  calculates  the  area  of  the  watershed,  and 
calculates  the  areal  rainfall  for  a  single  event.  The  second  program  (Program  B)  calculates  areal  rainfall  from  Equation 
9  for  a  series  of  events  after  first  calculating  the  weight  of  each  station.  Both  programs'  source  code  listings  are  in 
Appendixes  A  and  B  with  flow  charts.  A  detailed  description  for  preparing  the  input  data  for  the  programs  is  given  in 
Appendix  C.  Only  the  internal  functions  of  the  programs  are  explained  here. 

For  convenience  of  description,  three  types  of  points  are  defined  in  the  grid  system.  A  grid  point  occurs  wherever  an 
x  grid  line  intersects  a  y  grid  line,  point  A,  as  shown  in  figure  1.  A  boundary  point  occurs  wherever  the  watershed 
boundary  intersects  a  grid  line,  point  B,  figure  1.  A  break  point  occurs  wherever  the  watershed  boundary  has  an  abrupt 
change  in  direction,  point  C,  figure  1. 

In  Program  A,  after  the  data  have  been  read  in,  rainfall  at  grid  points,  boundary  points,  and  break  points  are 
calculated  from  Equation  1  and  results  are  stored  for  later  use. 

The  program  next  examines  each  square  mesh  in  turn  starting  from  the  origin  and  moving  up  the  first  column,  then 
back  to  the  bottom  of  the  next  column  to  the  right,  and  so  forth.  A  check  is  made  of  which  points  are  included  in 
each  mesh  and  which  of  these  points  are  inside  the  watershed  boundary. 

Starting  from  the  lower  left-hand  side  grid  point  of  a  square  mesh,  the  program  determines  whether  this  point  is 
located  inside  or  outside  the  watershed  boundary.  Only  the  inside  points  are  registered  for  storage.  Then  a  check  is 
made  for  the  existence  of  boundary  points  on  the  left-hand  sideline  of  the  mesh.  If  there  are  any,  they  will  also  be 
registered.  These  processes  are  repeated  for  the  next  grid  point  and  the  next  sideline  of  the  mesh  in  a  clockwise 
direction  until  all  four  grid  points  and  four  sides  of  square  mesh  have  been  examined,  if  all  the  points  are  located 
outside  the  watershed  boundary,  the  mesh  will  be  dropped  by  assigning  the  area  a  value  of  zero,  if  all  the  grid  points 
are  inside  the  watershed  boundary,  there  should  be  no  boundary  points  in  this  mesh  and  the  area  is  assigned  a  value  of 
one.  From  the  previous  point  rainfall  calculations,  the  average  volume  of  rainfall  for  this  mesh  can  be  quickly 
computed.  If  boundary  points  have  been  detected  in  the  mesh,  it  is  necessary  to  check  any  break  points  existing  in  the 
mesh.  If  break  points  are  included,  they  are  inserted  between  the  two  boundary  points  so  that  the  clockwise  order  of 
points  is  maintained. 
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By  using  all  the  registered  points  in  Equation  2,  the  area  of  a  boundary  mesh  can  be  obtained,  and  it  is  always  given 
as  a  fraction  of  one  since  the  area  of  a  complete  square  mesh  inside  the  watershed  was  assigned  as  one.  The  average 
rainfall  in  the  boundary  mesh  can  also  be  computed  from  the  information  available  in  the  previous  point  rainfall 
calculations.  After  every  mesh  in  the  grid  system  has  been  examined  and  calculated,  the  area  of  the  watershed  as  well 
as  the  areal  rainfall  can  be  computed  by  summing  the  individual  meshes. 

Most  parts  of  the  second  program  are  similar  to  that  of  Program  A  except  that  Program  B  does  not  calculate  point 
rainfall  and  the  average  rainfall  of  a  mesh.  Instead,  the  geometric  weights  (coefficient  part  of  Equation  8)of  each  mesh 
are  distributed  and  accumulated  into  the  corresponding  gages  to  compile  weights  given  in  Equation  9.  Once  a  mesh  is 
selected,  Program  B  counts  the  number  of  apexes,  mp,  and  calculates  the  area,  AA.  After  picking  up  a  point  in  the 
mesh,  the  program  calculates  the  fraction  of  the  reciprocal  squared  distance  of  nearby  n  gages,  multiplies  the  result  by 
AAp/mp,  and  accumulates  the  results  for  the  proper  gages.  By  repeating  the  same  procedure  for  all  points  and  meshes, 
the  value  accumulated  for  each  gage  station  divided  by  the  watershed  area  will  give  the  weight  of  each  station.  The 
program  then  will  read  a  rainfall  event  and  compute  areal  rainfall  from  Equation  9.  Since  the  last  procedure  can  be 
repeated,  this  program  can  calculate  as  many  rainfall  events  as  one  desires. 

DETERMINATION  OF  MESH  SIZE 
AND  NUMBER  OF  GAGES  USED  IN  CALCULATION 

In  determining  the  mesh  size  and  the  number  of  gages  to  be  used  in  calculating  rainfall  at  a  grid  point,  accuracy  and 
the  amount  of  computation  needed  should  be  the  main  concern.  Accuracy  of  computed  areal  rainfall  generally 
increases  with  an  increase  in  the  number  of  grid  points,  but  computer  time  increases  accordingly.  Therefore,  selection 
of  the  size  of  the  grid  system  is  a  compromise  between  the  cost  and  accuracy  of  computation.  More  points  will  give  a 
better  picture  of  the  rainfall  pattern  and  more  accurate  isohyets,  but  a  calculation  of  areal  rainfall  alone  (not  from  an 
isohyetal  map,  but  from  integrating  grid  point  rainfall)  would  need  less  points  than  the  isohyetal  map  without 
reducing  accuracy  significantly. 

The  number  of  gages  used  to  calculate  rainfall  at  a  grid  point  could  vary  from  two  gages,  the  minimum  number 
required  for  applying  Equation  1,  to  all  the  gages  in  the  network.  Since  the  summation  of  each  fraction  of  the 
reciprocal  distance  square  at  all  gages  always  totals  one  regardless  of  the  number  of  gages  used  in  the  calculation, 
increasing  the  number  of  gages  reduces  the  weight  of  nearby  gages  and  the  amount  of  rainfall  computed  tends  to 
become  more  or  less  an  average  value.  The  effect  of  the  local  distribution  pattern  of  the  rainfall  will  be  diminished. 
Also,  increasing  the  number  of  rain  gages  increases  computation  time.  Therefore,  only  a  few  of  the  available  gages 
should  be  used  at  each  grid  point  to  preserve  the  local  rainfall  distribution  pattern  as  well  as  to  minimize  the 
calculation.  Computer  resource  units,  an  index  of  the  amount  of  computation  involved  in  each  case,  are  listed  in  table 
1  to  give  examples  of  cost  which  will  depend  upon  suppliers  charge  per  CRU.  This  cost  does  not  include  computer 
storage  charges,  terminal  connect  time,  or  input/output  character  charges. 

The  storm  of  July  4-5,  1969,  in  the  Little  Mill  Creek  Watershed,  was  used  to  test  the  effects  of  mesh  size  and  the 
number  of  gages  in  calculating  grid  point  rainfall.  Mesh  sizes  of  2,000  ft.,  4,000  ft.,  6,000  ft.,  8*000  ft.,  and  10,000  ft. 


Table  1.  -  Computer  resources  units  used  in  calculation 


Number  of 
gages 

Number  of  grid  points 

143 

42 

25 

16 

16 

2  

9.97 

6.35 

5.60 

5.31 

5.35 

3  

11.13 

7.03 

6.09 

5.74 

5.78 

12.24 

7.58 

6.57 

6.19 

6.25 

5  

13.37 

8.20 

6.95 

6.59 

6.68 

6  

14.46 

8.73 

7.40 

7.01 

7.11 

8  

16.72 

9.85 

8.34 

7.82 

7.93 

10  

18.64 

10.91 

9.14 

8.55 

8.65 

12  

20.76 

11.93 

9.95 

9.28 

9.39 

5 


were  used  which  produced  143,  42,  25,  16,  and  16  grid  points  around  and  inside  the  watershed.  Also  the  number  of 
gages  used  in  computing  grid  point  rainfall  varied  from  a  minimum  of  two  to  all  12  gages.  The  location  of  rain  gages 
with  2,000  ft.  grid  system  on  the  Little  Mill  Creek  Watershed  is  given  in  figure  2. 

Areal  rainfall  calculated  from  the  combinations  of  different  arrangements  are  presented  in  table  2  Areal  rainfall 
calculated  from  the  isohyetal  map,  Thiessen  weight,  and  arthmetic  mean  were  5.52  inches,  5  49  inches  and  5  42 
inches  respectively  As  the  number  of  gages  used  at  each  grid  point  increased,  the  areal  rainfall,  regardless  of  the  mesh 
size,  decreased  and  became  closer  to  the  arithmetic  mean  as  expected. 


hie  2.  -  Areal  storm  rainfall  calculated  from  different  combinations  of  mesh 
and  number  of  gages  used,  Little  Mill  Creek,  July  4-5,  1969 


Number  of 

Areal  rainfall1  for 

grid  system 

spacing  of 

gages 

2,000  feet 

4,000  feet 

6,000  feet 

8,000  feet 

10,000  feet 

Inches 

Inches 

Inches 

Inches 

Inches 

2  

5.50 

5.51 

5.55 

5.52 

5.62 

3  

5.49 

5.49 

5.53 

5.51 

5.60 

4  

5.48 

5.50 

5.51 

5.51 

5.58' 

5  

5.47 

5.48 

5.50 

5.48 

5.57 

6  

5.46 

5.46 

5.48 

5.47 

5.55 

8  

5.44 

5.44 

5.46 

5.44 

5.53 

10  

5.43 

5.43 

5.44 

5.42 

5.51 

12  

5.41 

5.41 

5.42 

5.41 

5.49 

1  Rainfall  from  isohyetal  map  =  5.52  inches. 
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If  all  gages  had  equal  weight,  i.e.  arithmetic  mean,  the  weight  of  each  gage  would  be  0.083  (It  12).  In  general,  gage 
weights  tended  toward  the  arithmetic  mean  value  as  the  number  of  gages  used  at  each  grid  point  increases  from  2  to 
12.  The  larger  the  deviation  from  the  mean  for  the  2  gage  case,  the  greater  the  magnitude  of  change  as  gage  numbers 
increase  to  12.  These  relationships  indicate  the  averaging  effect  of  using  many  gages  in  the  areal  rainfall  calculation.  If 
we  assume  that  the  areal  rainfall  calculated  from  the  isohyetal  map,  5.52  inches,  is  correct,  mesh  sizes  of  2,000  ft.  and 
4,000  ft.  had  better  solutions  when  two  gages  were  used  (table  2),  and  6,000  ft.  mesh  size  for  three  or  four  gages.  The 
8,000  ft.  and  10,000  ft.  mesh  size  could  not  be  judged  for  reasons  that  will  be  mentioned  later.  Two  gages  were, 
therefore,  recommended  in  the  calculation  to  preserve  the  local  distribution  pattern  and  to  minimize  the  computer 
time. 

The  effect  of  mesh  size  in  table  2  is  not  as  clear  as  that  of  the  number  of  gages.  Mesh  sizes  of  2,000  ft.  and  4,000  ft. 
show  similar  results,  and  mesh  sizes  of  greater  than  6,000  ft.  do  not  show  a  clear  cut  pattern. 

When  a  coarse  grid  system  is  imposed  on  a  watershed,  such  as  8,000  ft.  and  10,000  ft.  mesh  size,  only  a  few  grid 
points  may  be  contained  inside  the  watershed  boundary.  In  this  case,  the  location  of  the  grid  points  can  become 
important  in  the  solution.  A  slight  change  in  the  grid  system  can  result  in  a  significant  change  in  the  solution.  As  the 
number  of  grid  points  increases,  the  location  of  the  grid  system  becomes  less  important.  Experience  suggests  that  the 
following  rule  may  be  used  in  determining  the  appropriate  mesh  size:  the  number  of  grid  points  which  fall  within  the 
watershed  area  should  be  about  three  to  four  times  the  number  of  rain  gages  available.  The  rule,  however,  should  not 
be  applied  on  a  dense  rain  gage  network  as  given  later  in  the  examples  for  Wayne  County  where  less  grid  points  are  not 
only  permissible  but  desirable. 

EVALUATION  OF  RDS  METHOD 

Three  gaged  watersheds  and  one  additional  area  were  selected  to  illustrate  the  application  and  the  validity  of  the 
RDS  method.  The  first  watershed,  No.  166  of  the  North  Appalachian  Experimental  Watershed,  Coshocton,  Ohio  is 
shown  in  figure  1.  The  area  of  WS-166  is  79.2  acres.  The  locations  of  five  rain  gages  are  shown  on  the  figure.  A  sixth 
gage  (113),  which  is  located  at  (4.70,  -6.60),  is  not  shown  on  the  figure. 

The  second  watershed  is  the  Little  Mill  Creek  Watershed  No.  97  of  Coshocton,  Ohio,  with  an  area  of  4,580  acres. 
The  boundaries  of  the  watershed  and  rain  gages  are  shown  in  figure  2. 

The  third  watershed,  No.  11,  shown  in  figure  3  is  part  of  the  Walnut  Gulch  watershed  in  Arizona.  It  was  used  by 
Diskin  and  Davis  (5)  in  their  study.  The  area  of  the  watershed  is  2,035  acres. 

A  dense  network  of  small,  tube-type  gages  has  been  organized  by  McGuinness  and  Harrold  (8)  over  an  area  of  1,500 
square  miles  in  northeast  Ohio  in  cooperation  with  1,700  volunteer  observers.  A  storm  measurement  in  Wayne 
County,  Ohio,  in  which  more  than  200  observations  were  reported,  were  used  in  this  report. 
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FIGURE  3.  —  Walnut  Gulch  watershed  No.  1 1  showing  rain  gage  locations  and  grid  system. 
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Areal  calculation 


The  area  of  a  watershed  was  determined  on  the  computer  by  summing  up  the  meshes  of  the  grid  system  which  was 
assigned  to  the  watershed.  Giving  the  area  of  a  square  mesh  a  value  of  one,  the  areas  of  the  watersheds  were  carefully 
counted  manually.  The  results  were  compared  with  computer  calculation  given  in  table  3.  Areas  of  the  watersheds  in 
acres  are  also  shown.  If  the  manual  calculation  is  accepted  as  the  true  area,  the  error  in  all  three  cases  is  less  than  1 
percent.  The  computer  calculation  always  gives  a  slightly  smaller  value  than  the  manual  method  since  the  convex 
boundary  lines  are  represented  by  a  chord.  The  accuracy  of  the  computer  calculated  area  can  be  increased  by 
increasing  the  number  of  break  points. 


Table  3.  —  Comparison  of  manual  and  computer  calculated  area 
[Unit  in  mesh  size] 


Watershed  name 

Area 

Manual 

Computer 

Error 

Acres 

Percent 

WS-166,  Coshocton  

79.2 

90.43 

90.23 

0.22 

4,580.0 

50.94 

51.18 

.47 

Walnut  Gulch,  11  

2,035.0 

80.85 

80.20 

.80 

Interpolation  of  missing  rainfall  records 

One  simple  application  of  the  RDS  method  is  interpolation  of  missing  rainfall  data  by  Equation  1.  A  commonly 
used  method  for  interpolation  is  that  of  the  U.S.  Weather  Bureau(9) .  Three  rain  gages  are  selected  as  evenly  spaced  as 
possible  around  the  rain  gage  where  the  record  was  missing.  The  normal  annual  rainfall  of  the  three  gages  is  compared 
to  that  of  the  data-missing  gage,  if  all  differences  are  less  than  ten  percent,  simple  arithmetic  mean  is  recommended. 
Otherwise,  the  normal-ratio  method  is  used  for  the  interpolation. 

Records  from  rain  gages  100,  103,  107,  Y101,  and  Y102  of  WS-166  (fig.  1)  were  used  in  an  interpolation  test.  It 
was  assumed  that  data  from  gage  100  were  missing.  Different  methods  were  used  to  interpolate  the  missing  data. 
Results  as  compared  with  measured  rainfall  are  given  in  table  4.  The  arithmetic  mean  and  normal-ratio  method  were 
calculated  from  rain  gages  103,  Y101  and  Y102  which  satisfies  the  Weather  Bureau  criteria.  All  gages  were  used  for 
interpolation  via  the  RDS  method.  All  methods  gave  satisfactory  results.  The  RDS  method  and  the  arithmetic  mean 
with  even  gage  distribution,  however,  show  a  closer  match  to  the  measured  data. 

The  normal-ratio  method  implies  that  the  storm  distribution  on  an  area  follows  the  distribution  pattern  of  normal 
annual  precipitation.  The  arithmetic  mean  method,  although  it  considers  each  event  individually,  assumes  that  all 
gages  have  the  same  weight  which  is  true  only  when  they  are  exactly  evenly  spaced.  The  RDS  method  not  only 
considers  individual  events  separately,  but  also  gives  more  weight  to  the  gage  that  is  closest  to  the  data  missing  station. 


Table  4.  —  Comparison  of  missing  rainfall  records  calculated  by  different  methods 


Date 

Measured 
rainfall 

RDS 

Arithmetic 
Even 
spacing1 

mean 
Uneven 
spacing2 

Normal 
ratio 

Inches 

Inches 

Inches 

Inches 

Inches 

April,  1967  

3.09 

2.96 

2.97 

2.87 

3.07 

July,  1967   

6.38 

6.55 

6.58 

6.52 

6.81 

March  28,  1967   

1.14 

1.06 

1.06 

1.05 

1.09 

June  8, 1967   

.17 

.17 

.17 

.18 

.18 

July  2,  1967   

1.36 

1.40 

1.42 

1.46 

1.47 

October  5,  1967   

.55 

.58 

.58 

.60 

.60 

1  Calculated  from  rain  gages  103,  Y101  and  Y102  (see  fig.  3). 
2 Calculated  from  rain  gages  103,  Y101  and  107  (see  fig.  3). 
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If  gage  107  is  substituted  for  gage  Y102  to  introduce  an  uneven  spacing,  the  results  become  worse  as  shown  in  table 
4.  This  indicates  that  the  RDS  method  gives  better  results  when  the  rain  gages  are  unevenly  spaced. 

Areal  rainfall  calculation  and  isohyetal  map 

The  storm  of  July  4-5,  1969  on  the  Little  Mill  Creek  watershed  was  used  for  the  areal  rainfall  calculation  and 
isohyetal  map  drawing.  The  isohyetal  map  drawn  by  conventional  methods  (manual  linear  interpolation)  is  given  in 
figure  4  as  solid  lines.  The  number  attached  to  each  gage  is  the  measured  rainfall.  Gage  MC  4  located  at  (7.08,  15.96) 
caught  6.61  inches  of  rainfall.  This  gage  was  not  shown  on  the  map  but  it  was  considered  in  drawing  the  isohyets. 

From  grid  point  rainfall,  the  isohyetal  map  was  drawn  manually  as  shown  in  figure  4  as  dotted  lines.  Although  the 
isohyetal  maps  drawn  by  the  two  different  methods  differ  in  detail,  the  general  pattern  of  both  maps  is  similar.  The 
areal  rainfall  calculated  by  the  RDS  method,  isohyetal  map,  Thiessen  weight  and  arithmetic  mean  are  given  in  the  first 
row  of  table  5.  All  methods  with  the  exception  of  the  arithmetic  mean  gave  about  the  same  result. 


  MANUAL 

FIGURE  4.  -  Isohyetal  map  for  storm  of  July  4-5,  1969,  Little 
Mill  Creek,  Coshocton,  Ohio. 


Table  5.  -  Areal  rainfall  calculated  from  different  methods 


Watershed 

Rainfall  period 

RDS 

Isohyet 

Thiessen 

Arithmetic 
mean 

Inches 

Inches 

Inches 

Inches 

Little  Mill  Creek 

July  4-5,  1969 

5.50 

5.52 

5.49 

5.42 

Walnut  Gulch  

July  28,  1966 

1.53 

1.53 

1.52 

1.63 

WS-166,  Coshocton  .  .  . 

1967 

34.53 

34.53 

34.48 

34.45 

April,  1967 

2.95 

2.94 

2.94 

2.95 

July,  1967 

6.41 

6.42 

6.40 

6.52 

March  28,  1967 

1.09 

1.09 

1.09 

1.07 

June  8,  1967 

.18 

.18 

.18 

.18 

July  2,  1967 

1.39 

1.41 

1.39 

1.43 

Oct.  5,  1967 

.59 

.59 

.59 

.59 
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Table  6.  —  Gage  weights  from  Thiessen  and  RDS  method,  Little  Mill  Creek 


Gage 

number 

Method 

MC4 

MC6 

27 

39 

54 

55 

56 

91 

Y101 

Y102 

Y103 

100 

Thiessen 

0 

0.01350 

0.23397 

0.17203 

0.07926 

0.06693 

0.12564 

0.23984 

0.02994 

0 

0.01350 

0.02539 

RDS 

.01350 

.02920 

.21392 

.16252 

.09921 

.06829 

.10168 

.22257 

.04243 

.00126 

.02258 

.02285 

The  gage  weights  for  each  rain  gage  calculated  from  the  Thiessen  and  the  RDS  method  are  given  in  table  6.  Gage 
MC4  and  Y102  have  no  effect  on  the  areal  rainfall  calculation  if  the  Thiessen  method  is  used  but,  in  the  RDS  method, 
both  gages  have  a  slight  weight. 

Diskin  and  Davis  (5)  used  a  computer  to  linearly  interpolate  the  measured  rainfall  at  two  rain  gages  and  drew  an 
isohyetal  map  for  the  storm  of  July  28,  1966,  on  the  Walnut  Gulch  watershed  (fig.  5).  Data  from  the  same  storm  were 
used  to  test  the  RDS  method.  The  isohyetal  map  based  on  the  grid  point  rainfall  is  shown  in  figure  6.  In  general,  the 
isohyetal  maps  drawn  by  the  two  different  methods  show  fair  agreement.  The  areal  rainfall  on  the  watershed,  as 
shown  in  the  second  row  of  table  5,  indicates  that  the  RDS  method  and  isohyetal  method  gave  almost  identical 
solutions  whereas  the  other  methods  gave  slightly  different  values. 

The  rain  gage  network  in  the  northeast  Ohio  area  measured  several  storms  during  its  operation.  One  of  them 
occurred  on  May  15,  1968,  with  a  storm  center  located  in  Wayne  County.  Two  hundred  and  seventy-five  observers 
reported  measurements  with  a  maximum  of  4.5  inches  in  the  storm  center.  A  manually  drawn  isohyetal  map  for  the 
storm  is  shown  as  solid  lines  in  figure  7.  In  the  calculation  of  grid  point  rainfall,  some  modification  was  made  in  the 
computer  program  to  reduce  calculation  time.  Instead  of  calculating  distances  from  a  grid  point  to  all  of  the  rain 
gages,  a  small  square  region  of  plus  and  minus  one  grid  unit  in  both  the  x  and  y  direction  from  the  grid  point  was 
delimited  and  only  rain  gages  in  that  region  were  used  to  calculate  the  rainfall  at  that  grid  point.  When  the  number  of 
rain  gages  included  in  the  region  was  fewer  than  the  number  of  rain  gages  used  to  calculate  the  grid  point  rainfall,  the 
region  was  expanded  by  one  grid  unit  until  the  necessary  number  of  rain  gages  was  included  in  the  region.  Because  of 
the  less  accurate  measuring  devices  used  in  the  rainfall  measurement  and  the  high  density  of  the  rain  gages,  four 
nearest  gages  were  used,  instead  of  two,  to  calculate  the  rainfall  at  a  grid  point.  The  isohyetal  map  from  grid  point 
rainfall  is  also  shown  in  figure  7  as  dotted  lines.  Area-depth  curves,  considering  only  the  area  inside  Wayne  County  are 
shown  in  figure  8.  The  solid  line  is  the  result  from  the  manually  drawn  isohyetal  map,  and  the  dotted  line  is  the  result 
from  the  isohyetal  map  drawn  from  grid  point  rainfall. 

Most  of  the  points  calculated  by  the  grid  point  method  gave  an  excellent  match  with  that  of  the  manual  method. 

In  a  further  investigation  of  the  proposed  method,  annual,  monthly,  and  daily  rainfall  for  1967  in  the  79-acre 
watershed,  166  shown  in  figure  1,  at  Coshocton  were  used.  Areal  rainfall  was  calculated  from  a  manually  drawn 


FIGURE  5.  -  Isohyetal  map  for  storm  of  July  28,  1966,  Walnut  FIGURE  6.  -  Isohyetal  map  drawn  from  grid  point  rainfall, 

Gulch  watershed,  Arizona  (from  Diskin  and  Davis,  1970).  Walnut  Gulch  watershed,  July  28,  1966. 
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isohyetal  map,  and  by  Thiessen  weight,  arithmetic  mean,  and  the  RDS  methods.  Results  are  presented  in  table  5  for 
annual  rainfall  of  1967,  monthly  rainfall  for  April  and  July  1967,  and  daily  rainfall  for  the  dates  shown  in  the  table. 

The  areal  rainfall  calculated  by  the  RDS  method  provides  a  closer  solution  to  that  of  the  isohyetal  method  than  the 
other  methods. 


RDS  METHOD 
MANUAL 


FIGURE  7.  —  Isohyetal  map  for  storm  of  May  15,  1968,  Wayne  County,  Ohio. 
 1  1  1  1  1  1  


2  5  10  20  50  100  200  500 

AREA,  SQUARE  MILES 

FIGURE  8.  —  Area-depth  curves  for  May  15,  1968,  storm  in  Wayne  County,  Ohio. 
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DISCUSSION 


The  Reciprocal  Distance  Square  method  is  an  attempt  to  make  the  best  use  of  existing  methods  to  develop  a 
computer-oriented  areal  rainfall  calculation  method.  Some  advantages  and  deficiencies  of  the  proposed  method  will  be 
discussed  in  this  section. 

In  this  study,  rainfall  at  a  grid  point  is  interpolated  from  the  measurement  of  nearby  gages  rather  than  calculated 
from  a  surface  trend  equation.  In  preliminary  trials,  the  bivariate  normal  equation,  suggested  by  Court  (2),  was  used 
to  calculate  areal  rainfall  as  well  as  grid  point  rainfall.  Although  the  method  provided  satisfactory  results  in  area-depth 
calculations,  it  was  unable  to  reproduce  a  satisfactory  isohyetal  map.  The  main  deficiency  of  trend-surface  fitting  is 
that  such  methods  fit  a  statistical  average  trend  of  the  rainfall  areal  distribution  pattern  at  the  expense  of  local 
distribution  tendency.  RDS  interpolation  only  considers  local  gages  so  that  the  area  distribution  patterns  were 
constructed  from  section  to  section  of  the  watershed  without  averaging  out  the  local  pattern.  Linear  interpolation  was 
also  tried,  but  the  results  were  unsatisfactory. 

Equation  1  is  a  rainfall  interpolation  equation  where  the  amount  of  rain  at  a  point  was  calculated  from  relative 
distances  from  the  point  to  the  nearby  gages.  It  is  also  possible  to  consider  other  factors  such  as  elevation  differences, 
slope,  and  orientation  of  topography  between  the  two  points.  It  is  not  the  purpose  of  this  study,  however,  to  develop 
such  an  interpolation  equation.  Rather,  the  purpose  is  to  show  that  if  the  rainfall  at  a  point  can  be  interpolated  from 
the  nearby  gages,  areal  rainfall  can  be  estimated  more  accurately  by  using  the  technique  given  in  this  report. 

The  satisfactory  results  obtained  in  the  current  study  indicate  that  distance  is  a  major  factor  in  a  small  areal  rainfall 
calculation.  In  applying  the  proposed  method,  therefore,  an  area  of  500  square  miles  is  suggested  as  a  limit.  Another 
restriction  should  be  on  mountainous  areas  where  the  elevation  of  each  point  might  have  a  significant  effect  on 
rainfall  calculation. 

Determination  of  mesh  size  and  the  number  of  gages  used  in  calcualtions  is  another  field  that  needs  further 
investigation.  The  selection  of  mesh  size  depends  on  several  factors,  such  as  the  size  of  the  area  or  the  watershed,  the 
number  of  rain  gages,  the  quality  of  rainfall  data,  the  type  of  computer  available,  and  the  purpose  of  the  study.  As 
shown  in  the  examples  given  in  the  previous  section,  too  coarse  a  grid  can  provide  an  unstable  solution.  A  finer  mesh 
avoids  such  problems  but  at  the  cost  of  increasing  the  amount  of  computation.  The  criteria  given  in  previous  sections 
can  be  used  as  a  general  guide  in  the  determination  of  the  mesh  size. 

Two  gages  were  used  to  compute  grid  point  rainfall  in  this  study  with  the  exception  of  the  dense  rain  gage  network 
in  Wayne  County  where  four  gages  were  used.  Although  a  fixed  number  of  gages  were  used  in  every  case,  it  is  possible 
to  use  a  variable  number  of  gages  by  assigning  a  radius  of  influence  and  using  all  the  gages  within  that  radius  to 
compute  rainfall  at  the  grid  point.  Of  course,  the  size  of  radius  is  another  question  that  will  arise  if  this  method  is 
used.  Technically,  this  method  is  available. 

The  most  important  advantages  of  the  new  method  are  objectivity  and  consistancy.  The  data  needed  for  the 
computation  are  rainfall  amounts,  coordinates  of  rain  gages,  boundary  intersection  points,  and  break  points. 
Coordinates  of  points  can  be  read  easily  without  any  human  judgement  being  involved.  All  the  calculations,  regardless 
of  whoever  prepares  the  data,  should  result  in  the  same  answer. 

Another  advantage  of  the  new  method  is  the  capability  of  utilizing  more  realistic  rainfall  input  to  many 
deterministic  watershed  models.  For  instance,  in  the  USDA  Hydrograph  Laboratory  Model/70  (7),  the  watershed  is 
grouped  into  several  zones  according  to  the  land-capability  classes  of  the  watershed.  Excess  rainfall  from  a  higher 
order  zone  is  cascaded  to  the  lower  order  zone.  With  the  help  of  grid  point  rainfall  calculated  by  the  new  method, 
rainfall  for  each  zone  can  be  computed  individually.  Or  if  desired,  each  zone  can  be  divided  into  several  sub-zones,  and 
the  rainfall  can  be  calculated  for  each  sub-zone  so  that  the  area  distribution  pattern  of  a  storm  can  be  retained  in  the 
watershed  model.  In  addition,  by  properly  assigning  the  geometric  factors  given  in  Equation  9,  it  is  possible  to 
calculate  the  gage  weight  for  each  zone.  Making  use  of  the  gage  weight  and  the  break  point  rainfall  data  in  rain  gages, 
the  break  point  rainfall  for  each  zone  can  be  composed.  The  input  data  to  the  model  can  preserve  the  important 
characteristics  of  rainfall— the  time  and  areal  distribution  pattern. 
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CONCLUSION 


A  new  approach  in  evaluating  precipitation  is  proposed.  Rainfall  at  a  point  is  assumed  to  be  proportional  to  the 
reciprocal  of  the  distance  squared  from  measuring  gages  to  the  point.  Methods  were  developed  to  calculate  arcal 
rainfall,  gage  weight,  grid  point  rainfall,  and  missing  data.  Two  computer  programs  were  written  to  c.irry  out  the 
calculations. 

The  new  RDS  method  was  applied  on  a  variety  of  data  where  the  area  ranges  from  a  few  acres  to  approximately  500 
square  miles  with  the  number  of  rain  gages  varying  from  six  to  several  hundred.  The  results  were  compared  with 
commonly  used  methods  including  the  isohyetal  map  method,  which  was  used  as  a  basis  for  comparison.  The  RDS 
method  was  found  to  be  close  to  that  of  the  isohyetal  method. 

The  RDS  method  with  the  help  of  the  computer  programs  not  only  reduced  the  manual  work  involved,  but  it  also 
provides  an  objective  means  of  calculating  areal  rainfall.  Once  the  gage  weights  for  each  rain  gage  have  been  calculated 
by  the  computer  program,  areal  rainfall  can  be  obtained  with  a  desk  calculator.  With  some  modification,  the  method 
can  estimate  a  partial  areal  rainfall  or  even  construct  an  areal  rainfall  histogram  to  express  the  temporal  and  spatial 
distribution  pattern  of  rainfall.  In  the  application  of  mathematical  watershed  models  on  the  hydrological  and 
environmental  problems,  these  data  should  provide  better  input  data  and  improve  the  solutions. 


LITERATURE  CITED 

(1)  CHIDLEY,  T.R.E.,  and  KEYS,  K.M. 

1970.  A  Rapid  Method  of  Computing  Area  Rainfall.  Jour,  of  Hydrol.  12(1):1  5-24. 

(2)  COURT,  A. 

1961.  Area-depth  Rainfall  Formulas.  Jour.  Geophys.  Res.  66(6):1823-1831. 

(3)  DISK1N,  M.H. 

1969.  Thicssen  Coefficients  by  Monte  Carlo  Procedure.  Jour,  of  Hydrol.  8(3):323-325. 

(4)   

1970.  On  the  Computer  Evaluation  ofThiessen  Weights.  Jour,  of  Hydrol.  ll(l):69-78. 

(5)    and  DAVIS,  P.R. 

1970.  A  Computer  Program  for  Interpolation  and  Plotting  of  Isohyetal  Points.  U.S.  Dept.  Agr.  ARS  41-177,  22  pp. 

(6)  HOLTAN,  H.N.  and  LOPEZ,  N.C. 

1971.  USDAHL-70  Model  of  Watershed  Hydrology,  U.S.  Dept.  Agr.  Tech.  Bui.  No.  1435,  84  pp. 

(7)  MANDEVILLE,  A.M.  and  RODDA,  J.C. 

1970.  A  Contribution  to  the  Objective  Assessment  of  Area  Rainfall  Amounts.  Jour,  of  Hydrol.  (N.Z.),  Wellington 
Symposium,  9(2):281-291. 

(8)  McGUINNESS,  J.L.  and  HARROLD,  L.L. 

1965.  Role  of  Storm  Surveys  in  Small  Watershed  Research.  Water  Resources  Res.  l(2):219-222. 

(9)  PAULHUS,  J.L.H.  and  KOHLER,  M.A. 

1952.  Interpolation  of  Missing  Precipitation  Records.  Mo.  Weather  Rev.  80:129-133. 

(10)  PENTLAND,  R.L.  and  CUTHBERT,  D.R. 

1971.  Operational  Hydrology  for  Ungaged  Stream  by  the  Grid  Square  Technique.  Water  Resources  Res.  7(2):283-291. 

(11)  SOLOMON,  S.I.,  DENOUVILLIEZ,  J.P.,  CHART,  E.J.,  WOOLEY,  J.A.,  and  CADOU,  C. 

1968.  The  Use  of  a  Square  Grid  System  for  Computer  Estimation  of  Precipitation,  Temperature,  and  Runoff.  Water 
Resources  Res.  4(5):919-929. 


13 


APPENDIX  A 
Computer  Program  A 
for  grid  point  rainfall  and  areal  rainfall  calculation 


START 


READ 
GENERAL 
INFOR- 
MATION, 


READ  WS 
BOUNDARY 
DATA 


READ 
GAGE 
LOCA- 
TION 


READ 
RAIN 
DATA 


CALCULATE  RAIN 
FALL  AT  ALL 
GRID  POINTS 


CALCULATE  RAIN- 
FALL AT  ALL 
BOUNDARY  POINTS 


CALCULATE  RAIN- 
FALL AT  ALL 
BREAK  POINTS 


SELECT 
A  MESH 


SET  AREA 
AND  RAINFALL 
TO  ZERO 


SET  AREA 
EQUAL  TO 
ONE 


CALCULATE 
PARTIAL 
AREA 


CALCULATE 

MEAN 

RAINFALL 

FOR 

AREA 

PRINT 
RESULTS 


STOP 


CALCULATE 
WATERSHED  AREA 

AND  AREA 
 RAINFALL 


FIGURE  9.  -  Flow  chart  for  calculating  watershed  area  and  rainfall,  Program  A. 
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C  THIS   PROGRAM  CALCULATES   THE   RAINFALL   AT  GRID  PO  IN  I  S  A\u. 

C         THE  AREA  RAINFALL 

C 

D I  MEMS  ION   YINT  (  22*  7 )*  XINT<22*  7  >*BREK(  1  00*  2  )  *  XG  (  30  )  ,  YG  (  30  )  * 
&RG<30)*HC22*  22)*  RYINTC22*  7 )*  RXINT (22*  7 )*  RBRK<  1 00 ) *  TEMPXC 20 ) * 
&TEMPY<20)*TEMPRC20)*KINDC20)*K5ID<20)*  AC  22* 22)*  RR( 22*  22) 

FILENAME  F1*F2*F3 

F1="WTRSH" 

F2="GAGE" 

F3="RA IN" 

DATA  NG AGE* N USED*  IMAX*  JMAX* I YMAX* JXM AX*  NBRK 
&/17*<l*  1  1*  13*  3*3*75/ 
C  READ  WATERSHED  DATA 

READCF1*300)  C  (YINT ( I  *  I Y )*  IY=1*7)*  I=1*IMAX) 
READCF 1  *  300 )  CCXINTC J* JX)*  JX=1*7)*  J=1*JMAX) 
READ(F1*301)    ( (BREK(N*  JB )*    JB=1*2)*  N=1*NBRK) 

300  FORMAT C 7F6 • 2 ) 

301  FORMAT ( 2F6 • 2 ) 

C         READ  RAINGAGE  LOCATION 

READ CF2* 302)    (XG(N)*   N= 1 *  NG AGE ) 

READ(F2*302)    <  YG (N ) *   N= 1  *  NG AGE ) 
C  READ  MEASURED  RAINFALL 

READ(F2*302)    (RG(N )*   N= 1 *NG AGE ) 

302  FORMAT ( ( 1 0F6 • 2 ) ) 
C         CLEAR  ARRAYS 

DO  12  I=1*IMAX 
DO  12  J«  1  *  JMAX 
R(  J* I )=0. 

12  CONTINUE 

DO    13   1=1* IMAX 
DO    13   IY= 1*1 YMAX 
RYINTC I* IY)=0. 

13  CONTINUE 

DO    14  J=1*JMAX 
DO    14  JX=1*JXMAX 
RXINT ( J* JX)=0. 

14  CONTINUE 

DO    17  N  =  1  *NBRK 
RBRK(N)=0. 
17  CONTINUE 
C         RAINFALL  AT  GRID  POINTS 
PRINT  210 

210  FORMAT C//2X* "CALCULATED  RAINFALL  AT  GRID  POINTS") 
DO    18    1=1* IMAX 

X=I 

TEMPXCL )=BREK ( JB* 1 ) 

211  FORMAT  C 1 H0 ) 

DO    18  J=1*JMAX 
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Y=  J 

RC  J* I )=RAIN (NGAGE* X, Y* XG*  YG* RG* MUSED) 
PRINT   212*  R(J*I) 
2  12  FORMATC 1H&*F5.2) 

18  CONTINUE 

C  RAINFALL   AT  BOUNDARY  POINTS 

PRINT  213 
DO    19  IMjIMAX 

IFCYINTCI*  1  ).LE.0.  )   GO   TO  19 
X=  I 

YINTCI* 1)=YINT(I* 1 )+l. 
I 1= Y INT ( I *  1  ) 
DO    19  IY=2*II 
Y=YINT( I» IY) 

RYINTC I* I Y)  =  RAIN (NG AGE* X*  Y»  XG*  YG»  RG>NUSED) 

19  CON  I INUE 

DO   22  J=1*JMAX 

IFCXINTC J* 1 ).LE.fl. )   GO   TO  22 
Y=  J 

XINTC J* 1 )=XINT( J* 1  )  +  l  . 
JJ=XINTC J> 1 ) 
DO  22  JX=2*JJ 
X=XINT( J*  JX) 

RXINTC J, JX)  =  RAIN (NGAGE*  X>  Y# XG* YG*  RG*NUSED) 

22  CONTINUE 

C         RAINFALL  AT  BREAKING  POINTS 
DO   23  N=1*NBRK 
X=BREKCN* 1 ) 
Y=BREKCN#  2 ) 

RBRKCN )  =  RAIN (NGAGE* X*  Y*  XG>  YG>  RG*NUSED> 

23  CONTINUE 

C  SELECT  POINTS   WITHIN   THE  WATERSHED 

AREA=0. 
RARE=0. 
JB=  1 

IMAX1=IMAX- 1 
JMAX1=JMAX- 1 
DO    190  I=1*IMAX1 
X=I 

XP=I+1 

DO    190  J=1*JMAX1 

Y=J 

YP=J+1 

L=0 

IFCYINTC I* 1 ).EQ. 0. )   GO   TO  30 
I  Y=2 

INT=YINTC I* 1 ) 
10  IYP=IY+1 

IFCY.LE.YINTC I, IY) )   GO   TO  20 
IFCY. GE. YINTC I* IYP) )   GO  TO  11 
L=L+1 
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TEMPX(L)=X 
TEMPYCL )= Y 
TEMPHCL)=HC J* I ) 
KINDCL)=0 
GO   TO  20 
11  IY=IY+2 

IFCIY.GT. INT)   GO   TO  30 
GO   TO  10 

20  IFCYINTCI* IY).LT. Y)   GO   TO  21 
IFCYINTC I* IY).GT. YP)   GO   TO  30 
L=L+  1 

TEMPXCL)=X 
TEMPYCL)  =  Y INT C  I  #  IY) 
TEMPRCL)  =  RYINT  C I  * IY) 
KINDCL)=1 

21  IFC IY.GE. INT)   GO   TO  30 
IY=I Y+l 

GO  TO  20 

30  IFCXINTC J+l, 1 ).EO.0. )   GO   TO  50 

JX=2 

JNT=XINTC J+l,  1 ) 

31  JXP=JX+1 

IFCX.LE- XINTC J* \, JX) )  GO   TO  40 
IFCX.GE.XINTC J+l, JXP) )  GO  TO  32 
L=L+  1 

TEMPX(L)=X 
TEMPY(L)=YP 
TEMPRCL)=RC J+l*  I  ) 
KINDCL )*0 
GO   TO  40 

32  JX=JX+2 

IFC JX.GT. JNT)  GO   TO  50 
GO  TO  31 

40  IFCXINTCJ+1, JX).LT.X)   GO   TO  41 
IFCXINTC J+l, JX).GT. XP)   GO  TO  50 
L=L+1 

TEMPXCL)=XINTC J+l, JX) 
TEMPYCL)=YP 
TEMPRCL)=RXINTC J+l,  JX) 
KINDCL)=1 

41  IFC JX.GE. JNT)  GO   TO  50 
JX=JX+1 

GO   TO  40 

50  IFCYINTC 1  +  1, 1 )»EQ. 0. )   GO  TO  70 
IY=YINTCI+1, 1 ) 

INT=2 

51  I YP= I Y- 1 

IFCYP.GE.YINTCI+1, IY))   GO   TO  60 

IFCYP.LE.YINTCI+1> IYP))  GO  TO  52 
L=L+  1 
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TEMPXCL)=XP 
TEMPYCL  >  =  YP 
TEMPRCL )=RC J+ 1  *  1  +  1  ) 
KINDCL )=0 
GO   TO   6  0 
52   I Y= I Y-2 

IF( IY.LT. INT)   GO   TO   7  0 
GO   TO   5  1 

60   IFCYINTC 1  +  1* IY).GT. YP )   GO   TO  6  1 
IFCYINTCI+1, IY).LT.  Y)   GO   TO  70 
L=L+  1 

TEMPXCL)=XP 
TEMPYCL)  =  YIiMTCI  +  l,  IY) 
TEMPRCL  )=RYINT  (  I  +  1*  IY) 
KINDCL)= 1 
6  1    IFC IY.LE. INT)   GO   TO  70 
I  Y=  I Y-  1 
GO   TO  6  0 

70  IFCXINTC J*  1  ).EQ.0.  )   GO   TO  90 
JX=XINTC J, 1 ) 

JNT=2 

71  JXP=JX-1 

IFCXP- GE. XINTC J, JX) )   GO   TO  80 
IFCXP-LE.XINTC J, JXP) )   GO  TO  72 
L=L«-1 

TEMPX(L)=XP 
TEMPYCL)=Y 
TEMPRCL )  =  RC  J* I  +  1 ) 
KINDCL)=0 
GO  TO  80 

72  JX*JX-2 

IFC JX.LT. JNT)  GO  TO  90 
GO   TO   7  1 

80  IFCXINTC J* JX).GT.XP)   GO   TO  81 
IFCXINTC J* JX).LT.X)   GO   TO  90 
L=L+  1 

TEMPXCL )SXINTC  J# JX) 

TEMPYCL )=Y 

TEMPRCL )=RXINTC J* JX) 

KINDCL)=1 

81  IFC JX.LE.  JNT)  GO   TO  90 
JX=  JX-  1 

GO  TO  80 

90  IFCL.GT. 1 )   GO   TO  91 
AC  J* I )  =  0. 

GO   TO  190 

91  IFCKINDC 1 ).EQ. 1 )   GO   TO  96 
DO   92  K=2>L 

IFCKINDCK).EQ.0)   GO   TO  92 
LOC=K+l 
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GO   TO  93 

92  CONTINUE 
AC  J* I )=  1. 
GO   TO  110 

93  IFC JB.GT.NBRK)   GO   TO  100 

IFCBREKC JB* 1 )• LT. X. 0 R. BREKC JB* 1 ).GT.XP)   GO   TO  100 
IFCBREKC J8» P).LT. Y.OR. BKEKC JB*2).GT. YP)   GO   TO  100 
DO   95  K=LOC*L 
KK=LOC+L+l-K 
KM=KK- 1 

TEMPXCKK)= I EMPXCKM) 
TEMP Y  C  KK )  =  TEM  PY  C KM ) 
K I N  D ( KK )  =  K I N  D ( KM ) 
TEMPRCKK)=TEMPRC KM ) 

95  CONTINUE 

94  TEMPXCL0C)=8R£KC JB* 1 ) 
TEMPYCLOC )  =  BREK( JB*  2  > 
KINDCL0C)=3 
TEMPRCLOC )  =  RBRKC  JB) 
L=L+  1 

JB=JB+1 
LOC=LOC+l 
GO  TO  93 

96  IFC JB.GT.NBRK)  GO   TO  100 

IFCBREKC JB*  1 ).LT.  X. OR. BREKC JB*  1 ).GT.XP)  GO  TO  100 
IFCBREKC JB*2).LT. Y.OR. BREKC JB*2).GT.YP)  GO  TO  100 
L=L+  1 

TEMPXCL )  =  BREKC  JB* 1 ) 

TEMPYCL)  =  BREKC  JB*  2) 

KINDCL)=3 

TEMPRCL )  =  RBRKC  JB) 

JB=JB+1 

GO   TO  96 
C       AREA  CALCULATION 
100  AC J* I )=FOFAC TEMPX* TEMPY* L ) 
110  AREA= AREA+ABS C AC  J* I ) ) 
C       RAINFALL  CALCULATION 

RRC J* I )s=FOFRC  TEMPR*  L ) 

RARE*RARE+RRC J* I )  *AC  J* I ) 
190  CONTINUE 

AMOUNT= RARE/ ARE A 

PRINT  201*  AREA 
201   F0RMATC//2X* "WATERSHED  AREA  =   "*F8.3*"  UNITS") 

PRINT  203*  AMOUNT 
203  FORMAT C//2X* "AMOUNT  OF   RAINFALL  ON   THE  WATERSHED  =  "*F3.3* 
&"INCHES") 

STOP 

END 

FUNCTION  FOFRCR* L ) 
C       CALCULATE  MEAN  RAINFALL   IN  DELTA  A 
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DIMENSION  R<20) 
XL=L 

P=0. 

DO    10    1=1 ,L 
P=P+R< I ) 
10  CONTINUE 
FOFR=P/XL 
HE TURN 
END 

FUNCTION  FOFACX*Y>N> 
C  CALCULATE   AREA  OF  POLYGON 

DIMENSION  X<20)>YC20> 
A=0. 
NM1=N-  1 
DO    10   1= 1 >  MM  1 

A=A+(X< 1+ 1 >*YC I >-X< I )*YC 1+1 ) ) 

10  CONTINUE 

A=A+X( 1 >*Y<N )~XCN >*YC 1 ) 
F0FA=A/2. 

RETURN 
END 

FUNCTION   RA  IN  CNO  Sis  XJ#  YI  *  X*  Y*  R*NUSED> 
C         CALCULATE  POINT  RAINFALL 

DIMENS ION   XC300)*  YC300)>  RC300)*  D<  300  )  *  Z  (  300  ) 

DEN=0. 

UP=0. 

DO    10  K=l*NOST 
DX=CX(K)-XJ)*(X(K)-XJ) 
DY=<YCK)-YI )*<Y<K)-YI ) 
DCK>=DX+DY 
Z<K)=»R<K) 
10  CONTINUE 

DO   20  K=l* MUSED 
KK=K+1 

DO   20  L-KK* NOST 

IFCD(K).LE.DCL) )   GO  TO  20 

TEMP=DCK) 

STOR=Z(K> 

D(K)=D(L) 

ZCK)=ZCL) 

D  (  L  >  =  T  EM  P 

Z(L)=STOR 
20  CONTINUE 

DO   30  K=1,NUSED 

UP=UP+Z(K)/DCK) 

DEN  =  DEN  +  1  •  /  D  (  K  ) 
30  CONTINUE 

RAIN=UP/DEN 

RETURN 

END 
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APPENDIX  B 
Computer  Program  B 
for  gage  weight  and  areal  rainfall  calculation 


READ 
GENERAL 

INFOR- 
MATION; 


READ 
BOUNDARY 
DATA 


READ  GAGE 
NAME  AND, 
LOCATE 


SELECT 
A  MESH 


SELECT  POINTS 
THAT  ARE  IN  SIDE 
THE  WATERSHED 


PICK  UP 
A  POINT 


CALCULATE  DIS- 
TANCE TO  ALL 
RAIN  GAGES 


SELECT  N  NEAR- 
EST RAIN  GAGES 
AND  DISTRIBUTE 
WT.  EACH  GAGE 


no 


PRINT  WT, 
FOR  ALL 
GAGES 


READ 
A  RAIN 
EVENT 


CALCULATE 
AREA  RAINFALL 


STOP 


FIGURE  10.  —  Flow  chart  for  calculating  station  weight  and  areal  rainfall  for  a  series  of  events,  Program  B. 
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C  THIS   PROGRAM   CALCULATES  GAGE  WEIGHT   AMD  AREA  RAINFALL 

C 

DIMENSION  YINT(22*7>*XINTC22*7>*BREKC100*2>*XGC30)*YG<30), 

*ftGC30>*MEFTC30)*DDC30>*WTGC30)*GN<30>*TEMPX<P0>*TEMPYCP0)* 
&KINDC20)* AC22*22> 

FILENAME  F1*F2*F3 

F  1  =  "WTRSH" 

F2  =  "6  AGE'* 

F3="RAIN" 

DATA  NG AGE*  N USED* IMAX*  JMAX* I YMAX*  JXMAX*  NBRK 
&/17*4* 1 1*  13*3*3*75/ 
C  READ  WATERSHED  DATA 

READCF1*300>  C < YINTC I* I Y)*  IY=1*7)*  I=1*IMAX> 
READ<F1*300)  CCXINTC J* JX>*  JX=1*7>*  J=1*JMAX) 
READ(F1*301)    C  (BHE'KCN*  «JB)*    JB=1*P>*  N=1*NBRK> 

300  FORMAT C 7F6 • 2 ) 

301  FORMAT C2F6. 2) 
DO    13   1=1* IMAX 

IFCYINTCI* 1 >*LE.0. )   GO   TO  13 
YINTCI* 1 )=YINT(I* 1 

13  CONTINUE 

DO    14  J=1*JMAX 

IF(XINTC J* 1 ).LE.0. )   GO   TO  14 
XINTC J* 1 >=XINT( J* 1 

14  CONTINUE 

C  READ  RAING AGE  NAME  AND  LOCATION 

READ(F2*302)  CXG(N>*  N=1*NGAGE> 
READ(F2*302)  ( YGCN )*  N=1*NGAGE) 
READCFP*303)    CGNCN )*   N= 1 *NGAGE ) 

302  FORMATC C 10F6.2) ) 

303  FO RM AT  < ( 1 5  A4  > ) 

C  SELECT   POINTS  WITHIN   THE  WATERSHED 

DO    12  N= 1  *  NG AGE 
WTGCN )=0. 
12  CONTINUE 
AREA=0. 
JB=1 

IMAX1=IMAX- 1 

JMAX1=JMAX-1 

DO    190  I=1*IMAX1 

X=I 

XP=I+1 

DO    190  J=1*JMAX1 

Y=  J 

YP=J+1 

L=0 

IFCYINTCI* 1 ).EQ.0.  )   GO   TO  30 
IY=2 
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INT-YINTC I*  1 ) 

10  IYP=IY+1 

IFCY.LE. YINTC I, IY) )   GO   TO  20 
IFCY-GE. YINTC l» IYP))   GO   TO  11 
L=L+  1 

TEMPXCL)=X 
TEMPYCL)=Y 
KINDCL)=0 
GO  TO  20 

11  IY=IY+2 

IFC IY.GT. INT)   GO   TO  30 
GO   TO  10 

20  IFCYINTC I, IY).LT. Y)   GO   TO  21 
IFCYINTC I* IY)«GT. YP)   GO   TO  30 
L=L+  1 

TEMPXCL)=X 

TEMPYCL )=YINTCI* IY) 

KINDCL)=1 

21  IFC IY.GE. INT)   GO   TO  30 
IY=IY+1 

GO  TO  20 

30  IFCXINTC J+ 1* 1 ).EQ. 0. )   GO   TO  50 
JX=2 

JNT=XINTC J+l* 1 ) 

31  JXP=JX+1 

IFCX.LE.XINTC J+l* JX) >   GO   TO  40 
IFCX.GE-XINTC J+l# JXP) )   GO   TO  32 
L«L+1 

TEMPX(L)=X 
TEMPYCL )=YP 
KINDCL)=0 
GO   TO  40 

32  JX=JX+2 

IFC JX.GT. JNT)   GO   TO  50 
GO  TO  31 

40  IFCXINTC  J+U  JX).LT.X)   GO  TO  41 
IFCXINTC J+l* JX).GT.XP)  GO   TO  50 
L»L+1 

TEMPXCL)  =  XINTC  J+U  JX) 

TEMPYCL)=YP 

KINDCL)*1 

41  IFC JX.GE. JNT)  GO   TO  50 
JX=JX+1 

GO   TO  40 

50  IFCYINTCI  +  U  1  ).EQ.0.  )  GO   TO  70 
IY=YINTCI+1# 1 ) 

IMT=2 

51  I YP= I Y- 1 

IFCYP.GE. YINTCI+ 1* IY))   GO  TO  60 
IF  CYP»LE»  YINTCI+1* IYP) )  GO   TO  52 
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L=L+  1 

TEMPXCL>=XP 
TEMPYCL )=YP 
KINDCL)=0 
GO   TO   6  0 
52  IY=IY-2 

IFC IY.LT. INT)   GO  TO  70 
GO   TO   5  1 

60  IFCYINTC I  +  l* IY).GT.  YP)   GO  TO  61 
IFCYINTCI+l* IY>.LT. Y>   GO   TO  70 
L=L+  1 

TEMPXCL)=XP 

T EM P Y CL)=YINTCI+1#  IY) 

KINDCD* 1 

61  IFC IY.LE. INT)   GO   TO  70 
IY=IY-1 

GO   TO  6  0 

70  IFCXINTC J* 1 ).EQ.  0.  )   GO   TO  90 
JX=XINT( J, 1 ) 

JNT=2 

71  JXP=JX-1 

IFCXP-GE.XINTC J,  JX) )  GO   TO  80 
IFCXP.LE.XINTC J* JXP) )   GO   TO  72 
L=L+1 

TEMPXCL)=XP 
TEMPYCL )=Y 
KIND(L>=0 
GO  TO  80 

72  JX=JX-2 

IFC JX.LT. JNT)   GO  TO  90 
GO  TO  71 

80  IFCXINTC J, JX).GT.XP)   GO   TO  81 
IFCXINTC J, JX).LT»X)   GO   TO  90 
L=L+1 

TEMPXCL)=XINTC J# JX) 

TEMPYCL)=Y 

KINDCL)=1 

81  IFC JX.LE. JNT)  GO   TO  90 
JX=  JX- 1 

GO  TO  80 

90  IFCL.GT.l)  GO  TO  91 
AC  J* I )  =  0. 

GO   TO  190 

91  IFCKINDC 1 ).EQ. 1 )   GO   TO  96 
DO  92  K=2*L 

IFCKINDCK).EQ.0)   GO  TO  92 

LOC=K+l 

GO  TO  93 

92  CONTINUE 
ACJ*I)=1. 
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GO   TO  110 

93  IF  <  JB«  GT • N  BRK )   GO   TO  100 

IF  C BREKC  JB* 1  )  •  LT«  X»  OR»  BREKC  JB*  1 ) • GT» XP )   GO    TO  100 
IFCBREKC  JB* 2).LT. Y.OR. BREKC  JB,2).GT.YP)   GO   TO  100 
DO   95  K=LOCL 
KK=LOC+L+l-K 
KM=KK- 1 

TEMPXCKK)=TEMPXCKM) 
TEMPYCKK)=TEMPYCKM> 
KINDC  KK)  =  KINDC  KM ) 

95  CONTINUE 

94  TEMPXCLOC)  =  BREKC JB,  1 ) 
TEMPYCLOC )  =  BREKC  JB#  2 ) 
KINDCLOC)=3 

L=L+1 
JB= JB+1 
LOC=LOC+l 
GO   TO  93 

96  IFC JB.GT.NBRK)   GO  TO  100 

IFCBREKC JB, 1 >.LT.X.OR. BREKC JB, 1 >.GT.XP>  GO  TO  100 
IFCBREKC JB*2).LT. Y. OR. BREKC JB, 2>.GT» YP)  GO  TO  100 
L=L+1 

TEMPXCL )= BREKC  JB>  1 ) 

TEMPYCL)=BREKC JB>2) 

KINDCL)=3 

JB= JB+1 

GO   TO  96 
C         CALCULATES  WATERSHED  AREA 
100  AC  J* I )=FOFACTEMPX#TEMPY*L> 
110  AREAS AREA+ABS  C  AC  J* I ) ) 
C         CALCULATES  GAGE  WEIGHT 

XL=L 

DO    120  K=1*L 

DO    121  N=1,NGAGE 

NEFT  CN )  =  0 

DX=CXGCN)-TEMPXCK) >*CXGCN>-TEMPXCK>  > 
DY=CYGCN)-TEMPYCK>  >*CYGCN  >-TEMPYCK> ) 
DDCN )=DX+DY 

121  CONTINUE 
KOUNT=l 

122  N=l 

123  IF C NEFT CN)»EQ»0)   GO  TO  124 
N»N  +  1 

GO   TO  123 

124  NEFT  CN )«  1 
NN*N 

TD=DDCN) 
MM=N+1 

DO    125  M=MM*NGAGE 
IFCNEFTCM).GT.0)   GO   TO  125 
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IFCTD.LE.DDCM>>   GO   TO  125 

N EFT CNN >  =  0 

NN=M 

NEFT  CM  >  = 1 
TD=DD(M) 

125  CONTINUE 
KOUNT=KOUNT+l 

IFCKOUNT.LE.NUSED>   GO   TO  122 
DEN=0. 

DO    126  N= 1 >  NG AGE 
IFCNEFTCN>.EQ.0)   GO   TO  126 
DEN=DEN+1  ./DDCN ) 

126  CONTINUE 

DO    127  N=  1*  NGAGE 

IFCNEFTCN>.EQ.0)   GO   TO  127 

UP=1./DDCN) 

DIS=UP/DFN 

?T=DIS/XL 

WTGCN>=WTGCN>+PT*AC J, I ) 

127  CONTINUE 
120  CONTINUE 
190  CONTINUE 

DO    130  N= 1 >  NGAGE 
WTG  CN )=WTG  CN )/ AREA 
130  CONTINUE 
PRINT  214 

214  FORMAT C//2X* "CALCULATED  WEIGHT  FOR  EACH  RAIN  GAGE") 
PRINT  215*    CGNCN  >*  WTGCN >*  N= 1 *NG AGE ) 

215  FO  RM  AT  C  5  X*  A4 .»  F  1  0  •  6  ) 

C  READ  MEASURED  RAINFALL 

READCF3>302)    C  RG  CN  >*   N= 1 *NG AGE ) 
C         CALCULATES  AREA  RAINFALL 

RAIN=0. 

DO    140  N= 1 > NGAGE 
RAIN=RAIN+RGCN)*WTGCN) 
140  CONTINUE 

PRINT  216*  RAIN 

216  F0RMATC//2X* "AMOUNT  OF  RAINFALL  =  M*F6»3> 
STOP 

END 

FUNCTION  FOFACX#Y*N> 
C         CALCULATES  AREA  OF  POLYGON 
DIMENSION  XC20)*YC20> 
A=0 

NM1=N- 1 

DO    10  I=1,NM1 

A=A+CXCI+1 >*YCI)-XCI)*YCI  +  1  )) 
10  CONTINUE 

A=A+XC 1 >*YCN )-XCN >*YC 1 ) 
F0FA=A/2. 
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APPENDIX  C 
Preparation  of  input  data 

When  the  programs  were  written,  some  restrictions  were  placed  on  the  input  data  to  keep  them  as  simple  as  possible 
to  avoid  confusion  and  to  save  labor.  It  is  strongly  recommended  that  anyone  planning  to  use  the  program  should  read 
this  section  thoroughly  and  carefully. 

The  following  materials  should  be  available: 

1.  A  watershed  map  that  indicates  the  watershed  boundary  and  available  rain  gages. 

2.  Graph  paper  of  a  size  to  cover  the  whole  watershed  map. 

3.  Rainfall  records  from  the  rain  gages. 

Before  tracing  the  watershed  map  on  the  graph  paper,  it  is  necessary  to  examine  the  meshes  that  contain  the 
watershed  boundary  so  that  the  following  situations  do  not  occur.  First,  there  should  not  be  two  or  more  separated 
watershed  areas  in  a  mesh  as  shown  in  figure  11-A  by  the  shaded  area.  Second,  if  possible,  limit  the  number  of 
boundary  points  in  the  mesh  to  no  more  than  two  to  avoid  forming  separted  outside  areas  in  a  mesh  as  shown  in 
figure  11-B. 

Violations  of  the  above  rules  can  be  eliminated  by  shifting  and  rotating  the  graph  paper  on  the  map.  The  first  rule 
should  always  be  obeyed  but,  if  the  second  rule  must  be  violated,  no  break  points  should  be  assigned  to  this  mesh. 
The  error  in  calculating  watershed  area  caused  by  this  situation  is  negligible. 

After  the  watershed  boundary  and  gage  stations  have  been  traced  on  graph  paper,  a  Cartesian  coordinate  system, 
with  origins  at  (1,  1)  is  drawn  so  that  the  whole  watershed  area  lies  within  the  first  quadrant.  The  lowest  point  of  the 
watershed  boundary  and  the  furthest  point  on  the  left  hand  side  of  the  watershed  boundary  should  come  as  close  as 
possible  to  the  axes.  The  gage  stations,  however,  can  be  located  in  any  quadrant.  A  grid  system  covering  the  entire 
watershed  area  can  then  be  drawn  as  shown  in  figures  3,  4,  and  5. 

The  input  data  required  manually  from  the  user  at  a  terminal  for  both  programs  are  almost  identical,  and,  unless 
otherwise  specified,  are  assumed  applicable  to  both.  The  first  set  of  data  input  contains  general  information  and  is 
given  in  the  program  as 

DATA  NGAGE,  NUSED,  IMAX,  JMAX,  IYMAX,  JXMAX,  NBRK 

where 

NGAGE  is  the  total  number  of  rain  gages  available  for  use, 

NUSED  is  the  number  of  rain  gages  that  are  going  to  be  used  in  the  calculation.  This  number  must  not  exceed 
NGAGE  and  values  of  2  to  4  are  recommended, 

IMAX  is  the  maximum  unit  of  the  x  grid  line, 

JMAX  is  the  maximum  unit  of  the  y  grid  line, 

IYMAX  is  the  maximum  number  of  boundary  points  on  any  x  grid  line, 
JXMAX  is  the  maximum  number  of  boundary  points  on  any  y  grid  line,  and 
NBRK  is  the  total  number  of  break  points. 

From  a  file,  Fl,  the  program  READs  information  concerning  the  watershed  boundary  with  these  statements: 

READ  (Fl,300)  ((YINT  (I,IY),  IY=1,7),  1=1,  IMAX) 

READ  (F1,300)((XINT(J,JX),  JX=1,7),  J=1,JMAX) 

READ  (Fl,301)  ((BREK(NJB),  JB=1,2),  N=1,NBRK). 

YINT  (I,IY)  are  the  boundary  points  where  the  boundary  line  crosses  the  x  grid  lines.  If  the  boundary  point  coincides 
with  a  grid  point,  it  is  regarded  as  a  boundary  point.  IY  denotes  the  y-coordinate  of  the  boundary  points  except  that 
the  first  space  (IY=1)  is  reserved  to  store  the  total  number  of  the  boundary  points  on  the  x  grid  line.  I  represents  units 
of  the  x-axis  and  runs  from  1=1  to  I=IMAX.  Within  a  grid  line,  boundary  points  are  arranged  in  an  ascending  order  of 
the  y  values.  If  no  boundary  points  exist  in  a  grid  line,  zero  should  be  used  for  at  least  the  first  column  (IY=1). 
XINT(J  JX)  are  the  points  where  y  grid  lines  intercept  the  watershed  boundary  and  should  be  prepared  in  the  same 
manner  as  YINT(I,IY). 
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FIGURE  11.  —  Examples  of  undesirable  situations  in  the  grid  system.  FIGURE  12.  —  Rules  to  assign  the  order  of  break  points. 

Break  points  are  read  in  by  the  last  read  statement  of  BREK(NJB).  In  selecting  the  break  points,  the  second  rule 
mentioned  in  the  early  part  of  this  section  should  be  observed.  If  a  grid  line  is  tangent  to  the  watershed  boundary,  the 
tangential  point  is  considered  as  a  break  point.  The  arrangement  of  the  data  is  as  follows.  If  there  is  more  than  one 
break  point  in  a  mesh,  as  shown  in  figure  12,  the  break  points  should  be  arranged  numerically  in  a  clockwise 
direction.  In  figure  12,  the  solid  lines  indicate  the  inside  of  the  watershed  and  the  dotted  lines  indicate  the  outside  of 
the  watershed. 

The  break  points  in  each  mesh  are  listed  consecutively  starting  with  the  lowest  mesh  of  the  left  most  column  of  the 
grid  system.  After  reaching  the  uppermost  mesh  of  the  column,  change  to  the  next  column  and  start  from  the 
lowermost  mesh  again.  Each  break  point  will  be  registered  in  one  line  with  x-coordinate  value  first  and  the 
y-coordinate  value  next. 

The  third  part  of  the  input  data  is  information  on  the  location  of  all  rain  gages  read  from  the  F2  file  by 

READ  (F2,  302)  (XG(N),  N=1,NGAGE) 

READ  (F2,  302)  (YG(N),  N=1,NGAGE). 

XG(N)  is  the  x-coordinate  and  YG(N)  is  the  y-coordinate  of  the  rain  gage.  No  order  was  assigned  for  arranging  the 
data  but  the  order  of  listing  XG(N)  and  YG(N)  must  correspond.  Each  line  contains  10  data  points  or  less.  In 
addition,  Program  B  also  needs  the  name  or  number  of  the  gage  stations.  These  are  read  from  file  F2  by 

READ  (F2,303)  (GN(N),  N=1,NGAGE). 

An  alphanumeric  field  of  four  characters  was  allotted  to  each  gage,  and  the  gages  should  be  listed  in  the  same  order  as 
the  rain  gage  coordinates.  No  more  than  15  data  can  be  put  on  a  line. 
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The  last  part  of  the  input  data  are  depth  of  measured  rainfall  at  each  gage.  The  data  can  be  annual,  seasonal, 
monthly,  weekly,  or  daily  rainfall,  or  rainfall  for  a  special  storm. 

For  Program  A,  data  are  read  from  file  F2  with 

READ(F2,302)  (RG(N),N=1  ,NGAGE)- 

to  read  one  event  only.  In  Program  B,  data  are  read  from  file  F3  with 

READ  (F3,302,END=999)  (RG(N),  N=1,NGAGE) 

for  continuous  reading  of  several  events.  The  data,  of  course,  should  be  arranged  in  the  same  order  as  the  other  rain 
gage  data. 

It  should  be  pointed  out  again  that  computer  capacity  was  a  factor  when  the  programs  were  developed.  Some 
modifications  of  the  program  will  be  required  if  it  is  to  be  used  on  a  larger  computer  or  for  a  larger  size  watershed. 
The  main  changes  occur  at  DIMENSION,  READ  and  PRINT  statements.  The  teletype  used  limited  output  to  72 
characters  per  line.  Other  output  devices  may  allow  up  to  132  characters  per  line.  Modification  of  the  format  is 
recommended  to  make  best  use  of  these  larger  input-output  sizes.  With  some  knowledge  of  programming  or  with  the 
help  of  a  consultant,  this  can  be  done  without  too  much  difficulty.  No  identification  data  such  as  the  name  of  the 
watershed  or  the  date  of  storm  were  given  in  the  programs.  If  desired  they  can  be  added  easily  with  several  extra 
statements. 

A  few  words  on  the  output  might  be  helpful  although  most  of  it  is  self-explanatory.  In  Program  A,  rainfall  at  grid 
points  is  given  in  the  output.  Grid  point  rainfall  is  printed  out  in  matrix  form  so  that  the  results  can  be  used  directly 
in  drawing  an  isohyetal  map.  An  example  of  the  output  is  given  in  figure  13.  The  areal  rainfall  is  given  so  that  it  is  not 
necessary  to  measure  the  isohyetal  map  and  calculate  the  areal  rainfall. 

In  Program  B,  the  results  are  more  straight  forward.  First,  it  prints  out  the  weight  of  each  rain  gage,  and  then  it 
prints  out  the  areal  rainfall  of  each  event. 


CALCULATED  RAINFALL  AT  GRID  POINTS 

4.56  4.54  4.53  4.56  4.94  5.27  5.64  5.69  5.70  5.68  5.66  5.93  6.03 
4.55    4.50    4.49    4.55    4.73    5.53    5.65    5.73    5.73    5.70  5.66  5.94  6.07 

4.57  4.50  4.48  4.59  4.85  5.42  5.59  5.71  5.72  5.69  5.69  5.99  6.13 
4.64  4.65  4.74  4.76  5.10  5.14  5.43  5.60  5.69  5.68  5.98  6.07  6.23 
4.71  4.90  4.88  4.94  5.17  5.20  5.28  5.33  5.61  5.62  6.15  6.20  6.63 
4.87  4.94  4.96  4.96  5.14  5.28  5.38  5.39  5.88  6.11  6.23  6.65  6.63 
4.77  4.86  4.97  5.07  5.07  5.40  5.45  5.59  6.01  6.31  6.39  6.66  6.64 
4.77  4.77  5.22  5.28  5.21  5.34  5.42  5.78  6.33  6.57  6.55  6.66  6.64 
4.76  4.93  5.13  5.13  5.28  5.29  5.47  6.11  6.54  6.68  6.62  6.67  6.65 
4.99  5.01  5.07  5.41  5.37  5.32  5.48  6.25  6.51  6.61  6.58  6.67  6.65 
4.99  5.00  5.00  5.32  5.34  5.32  5.31  6.27  6.42  6.50  6.51  6.67  6.66 
AMOUNT  OF  RAINFALL  ON  THE  WATERSHED  =  5.50  INCHES 

PROGRAM  STOP  AT  1445 

USED     8.07  UNITS  (Computer  Resource  Units) 

Note:  Orientation  of  this  output  was  rotated  90  degrees  of  that  of  figure  4. 

FIGURE  1 3.  -  Computer  output  of  grid  point  rainfall,  Little  Mill  Creek,  July  4-5,  1969. 
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